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Abstract 

We propose a full-wave pseudo-analytical numerical electromagnetic (EM) algorithm to model subsurface induc¬ 
tion sensors, traversing planar-layered geological formations of arbitrary EM material anisotropy and loss, which are 
used, for example, in the exploration of hydrocarbon reserves. Unlike past pseudo-analytical planar-layered modeling 
algorithms that impose parallelism between the formation’s bed junctions however, our method involves judicious 
employment of Transformation Optics techniques to address challenges related to modeling relative slope (i.e., tilt¬ 
ing) between said junctions (including arbitrary azimuth orientation of each junction). The algorithm exhibits this 
flexibility, both with respect to loss and anisotropy in the formation layers as well as junction tilting, via employing 
special planar slabs that coat each “flattened” (i.e., originally tilted) planar interface, locally redirecting the incident 
wave within the coating slabs to cause wave fronts to interact with the flattened interfaces as if they were still tilted 
with a specific, user-defined orientation. Moreover, since the coating layers are homogeneous rather than exhibiting 
continuous material variation, a minimal number of these layers must be inserted and hence reduces added simulation 
time and computational expense. As said coating layers are not reflectionless however, they do induce artificial field 
scattering that corrupts legitimate field signatures due to the (effective) interface tilting. Numerical results, for two 
half-spaces separated by a tilted interface, quantify error trends versus effective interface tilting, material properties, 
transmitter/receiver spacing, sensor position, coating slab thickness, and transmitter and receiver orientation, helping 
understand the spurious scattering’s effect on reliable (effective) tilting this algorithm can model. Under the effective 
tilting constraints suggested by the results of said error study, we finally exhibit responses of sensors traversing three¬ 
layered media, where we vary the anisotropy, loss, and relative tilting of the formations and explore the sensitivity of 
the sensor’s complex-valued measurements to both the magnitude of effective relative interface tilting (polar rotation) 
as well as azimuthal orientation of the effectively tilted interfaces. 

Keywords: multi-layered media, deviated formations, geological unconformity, induction well-logging, borehole 
geophysics, geophysical exploration 


1. Introduction 

Long-standing and sustained interest has been directed towards the numerical evaluation of electromagnetic (EM) 
fields produced by sensors embedded in complex, layered-medium environments (T). In particular, within the con¬ 
text of geophysical exploration (of hydrocarbon reserves, for example), there exists great interest to computationally 
model the response of induction tools that can remotely sense the electrical and structural properties of complex 
geological formations (and consequently, their hydrocarbon productivity) mm. Indeed, high-fidelity, rapid, and 
geometry-robust computational forward-modeling aids fundamental understanding of how factors such as formation’s 
global inhomogeneity structure, conductive anisotropy in formation bed layers, induction tool geometry, exploration 
borehole geometry, and drilling fluid type (among other factors) affect the sensor’s responses. This knowledge informs 
both effective and robust geophysical parameter retrieval algorithms (inverse problem), as well as sound data interpre¬ 
tation techniques ESI. Developing forward-modeling algorithms which not only deliver rapid, accuracy-controllable 
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results, but also simulate the effects of a greater number of dominant, geophysical features without markedly increased 
computational burden, represents a high priority in subsurface geophysical exploration dum. 

In the interest of obtaining a good trade-off between the forward modeler’s solution speed while still satisfactorily 
modeling the EM behavior of the environment’s dominant geophysical features, a layered-medium approximation of 
the geophysical formation often proves very useful. Indeed cylindrical layering, planar layering, and a combination 
of the two (for example, to model the cylindrical exploratory borehole and invasion zone embedded within a stack 
of planar formation beds) are arguably three of the most widely used layering approximations in subsurface geo¬ 
physics mmrnmm, for both onshore and offshore geophysical exploration modeling Il23]j33ll . The prevalence of 
layered-medium approximations stems in large part, at least from a computational modeling standpoint, due to the 
typical availability of closed-form eigenfunction expansions to compute the EM field (34|[Ch. 2-3]. These full-wave 
techniques are quite attractive since they can robustly deliver rapid solutions with high, user-controlled accuracy under 
widely varying conditions with respect to anisotropy and loss in the formation’s layers, orientation and position of the 
electric or (equivalent) magnetic current-based sensors (viz., electric loop antennas), and source frequency USED. 
The robustness to physical parameters is highly desirable in geophysics applications since geological structures are 
known to exhibit a wide range of inhomogeneity profiles with respect to conductivity, anisotropy, and geometrical 
layering MHO). For example, with respect to formation conductivity properties, diverse geological structures can 
embody macro-scale conductive anisotropy in the induction frequency regime, such as (possibly deviated) sand-shale 
micro-laminate deposits, clean-sand micro-laminate deposits, and either natural or drilling-induced fractures. In the 
sub-2MHz frequency regime, the electrical conduction current transport characteristics of such structures indeed are 
often mathematically described by a uniaxial or biaxial conductivity tensor exhibiting directional electrical conduc¬ 
tivities whose value range can span in excess of four orders of magnitude ll2l 17] [12Tl . 

When employing planar and cylindrical layer approximations one almost always assumes that the interfaces are 
parallel, i.e. exhibit common central axes (say, along z ) in the case of cylindrical layers mm, or interfaces that are 
all parallel to a common plane in the case of planar layers 13 Ho). However, it may be more appropriate in many cases 
to admit layered media with material property variation along the direction(s) conventionally presumed homogeneous. 
For example, in cylindrically-layered medium problems involving deviated drilling, gravitational effects may induce 
a downward diffusion of the drilling fluid that leads to a cylindrical invasion zone angled relative to the cylindrical 
exploratory borehole E). Similarly, formations that locally (i.e., in the proximity of the EM sensor) appear as a 
“stack” of beds with tilted (sloped) planar interfaces can appear (for example) due to temporal discontinuities in 
the formation’s geological record. These temporal discontinuities in turn can manifest as commensurately abrupt 
spatial discontinuities, known as unconformities (especially, angular unconformities) DSH23 ; see Figure [Ta] below 
for a schematic illustration. Indeed, the effects of unconformities and other complex formation properties (such as 
fractures) have garnered increasing attention over the past ten years |[38l|40l , particularly in light of the relatively 
recent availability of induction sensor systems offering a rich diversity of measurement information with respect 
to sensor radiation frequency, transmitter and receiver orientation (“directional” diversity), and transmitter/receiver 
separation ED}®. 

A natural question arises as to which numerical technique is best suited to modeling these more complex geome¬ 
tries involving tilted layers. In principle, one could resort to brute-force techniques such as finite difference and finite 
element methods (3 HU 53 53 IS) • The potential for low-frequency instability (e.g., when modeling geophysical 
sensors operating down to the magnetotelluric frequency range (fraction of a Hertz)), high computational cost (unac¬ 
ceptable especially when many repeated forward solutions are required to solve the inverse problem), and accuracy 
limitations due to mesh truncation issues (say, via perfectly matched layers or other approaches B71148) ) associated 
with the lack of transverse symmetry in the tilted-layer domain ||49l , render these numerical methods less suitable for 
developing fast forward-modeler engines for tilted-layer problems. Another potential approach involves asymptotic 
solutions which traces the progress of incident rays and their specular reflections within subsurface formations ll50l . 
However, this approach is limited to “sufficiently high” frequencies and hence may often be unsuitable for modeling 
induction-regime sensors operating in zones where highly resistive and highly conductive (not to mention anisotropic) 
layers may coexist (ED. Yet a third possible approach, called the “Tilt Operator” method, which assumes lossless me¬ 
dia and negligible EM near-fields to avoid spurious exponential field growths (arising from violation of “primitive” 
causality [i.e., cause preceding effect], which is inherent in this method), is another possibility Ii52ll53) . Akin to the 
other mentioned high-frequency approach however l50l . the Tilt Operator method is not appropriate for our more 
general class of problems with respect to sensor and geological formation characteristics. 
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Turning our attention henceforth to tilted planar-layered media, we propose a pseudo-analytical method based on 
EM plane wave eigenfunction expansions that manifest mathematically as two-dimensional (2-D) Fourier integrals. 
This is in contrast to faster, but more restrictive (with respect to allowed media) 1-D Fourier-Bessel (“Sommerfeld”) 
and Fourier-Hankel integral transforms that express EM fields in planar-layered media as integral expansions of EM 
conical wave eigenfunctions lf34l lCh. 2]. Our choice rests upon robust error control capabilities and speed perfor¬ 
mance of the 2-D integral transform with respect to source radiation frequency, source distribution, and material prop¬ 
erties 0E3. The use of eigenfunction expansions for modeling EM behavior of non-parallel layers is enabled here by 
the use of Transformation Optics (T.O.) techniques ll47l [54l{58l to effectively replace the original problem (with tilted 
interfaces) by an equivalent problem with strictly parallel interfaces, where additional “interface-flattening” layers 
with anisotropic response are inserted into the geometry to mimic the effect of the original tilted geometry (c.f. Fig. 
[Tb|. We remark that exploitation of T.O. techniques to facilitate numerical field computation based on brute-force al¬ 
gorithms—rather than the eigenfunction-based spectral approach considered here—was recently proposed in If59li60il . 
We moreover emphasize that the 2-D Fourier integral is capable of modeling propagation and scattering behavior of 
these interface-flattening media, which possess azimuthal non-symmetric material tensors; the two stated 1-D integral 
transforms, which are inherently restricted to modeling azimuthal-symmetric media, by contrast lack this numerical 
modeling capability. 

Before proceeding, we note that the coating slabs are spatially homogeneous in the employed Cartesian coordinate 
system, and hence we require only one planar layer to represent each coating slab, and that too to represent each slab’s 
spatial material profile exactly. This homogeneity characteristic is important from a computational efficiency stand¬ 
point, as it means that for each “flattened” interface only two coating layers’ (source-independent) EM eigenfunctions, 
Fresnel reflection and transmission matrices, etc. need to be computed ||T9l . Second, as we will be fundamentally ap¬ 
proximating the transverse translation-variant geometry as a transverse translation-invariant one, modeling spurious 
scattering from the “apexes” and more complex intersection junctions of the tilted beds is out of the question. As we 
concern ourselves with subsurface geophysical media, which typically present inherent conductivity (typically on the 
order of at least 10 1 * 3 S/m to 2S/m J2j|), scattering from these intersections should typically be negligible so long as 
the sensors are not in the immediate neighborhood of said intersections (rarely the case, for the small tilting explored 
herein). Third, the interface-flattening slabs are not impedance-matched to the respective ambient medium layers into 
which they are inserted. Indeed, one of the objectives of this paper is to quantify the impedance mismatch of the 
artificial slabs, which we will find practically constrains (i.e., for a given desired computation accuracy) the amount 
of interface tilt that can be reliably modeled. 

This paper is organized as follows. In Section [2] we overview the 2-D plane wave expansion algorithm, derive 
the material blueprints for the planar “interface-flattening” coating slabs, and show how to systematically incorporate 
these into the computational model. Sections 3.1|3.2 show the error analysis to quantify how the accuracy of the 
results varies with effective interface tilting, material profile, transmitter/receiver spacing, sensor position, coating 
slab thickness, and complex-valued measurement component (both its real and imaginary parts). In Section [4] we 
apply the algorithm to predicting EM multi-component induction tool responses when the tool traverses (effectively) 
tilted formation beds for different interface tilt orientations as well as central bed anisotropic conductivity profiles. The 
formation anisotropies studied will span the full gamut: All the way from isotropic to (“Transverse-Isotropic”) non- 
deviated uniaxial, (“cross-bedded”) deviated uniaxial, and full biaxial anisotropy. We adopt the exp (-iu>t) convention, 
as well as assume all EM media are spatially non-dispersive, time-invariant, and are representable by diagonalizable 
anisotropic 3x3 material tensorsQ 


2. Formulation 

2.1. Background: Electromagnetic Plane Wave Eigenfunction Expansions 

In deriving the planar multi-layered medium eigenfunction expansion expressions, first assume a homogeneous 
formation whose dielectric (i.e., excluding conductivity), relative magnetic permeability, and electric conductivity 


1 Diagonalizability of the material tensors, which physically corresponds to a medium having well-defined electrical and magnetic responses for 

any direction of applied electric and magnetic field, is required for completeness of the plane wave basis. All naturally-occurring media, as well as 

the introduced interface-flattening slabs applied to geophysical media, are characterized by diagonalizable material tensors. 
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constitutive anisotropic material tensors write as e r , p r , and <f. Specifically, the assumed material tensors are those 
of the layer (i.e., in the anticipated multi-layered case), labeled M, within which the transmitter resides. Maxwell’s 
Equations in the frequency domain, upon impressing causative electric current Jf (r) and/or (equivalent) magnetic 
current Wl(r), yields the electric field vector wave equation (duality in Maxwell’s Equations yields the magnetic field 
vector wave equation) mmE 

M-) - V x p r l ■ V x -kl (e r + id-/cj) •, 3K (£) = ikorjoST - V x JT r l ■ M (2.1) 


Now define the three-dimensional spatial Fourier Transform (FT) pair for some generic vector field _£ (e.g., the 
magnetic field or current source vector) E): 


L(k) 


+oo 




dxdydz, -£(r) 



-oo 


d k x dk v dk z 


( 2 . 2 ) 


where r = (x,y, z) is the position vector and k = (k x ,k y ,k z ) is the wave vector. Expanding the left and right hand 
sides, of the second equation in Eqn. in their respective wave number domain 3-D integral representations and 
matching the Fourier-domain integrands on both sides, one can multiply the inverse of A (the FT of JL) to the left of 
both integrands. Admitting a single Hertzian/infinitesimal-point transmitter current source located at r' = (x',y', z!), 
and denoting the receiver location r, one can then procure the “direct” (i.e., homogeneous medium) radiated time- 
harmonic electric field £,/(r) |fl9l . Indeed, performing “analytically” (i.e., via contour integration and residue calculus 
techniques) the k z integral leads to the following expression: 


Too 
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ikM, nz Az 
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ikM.nz^Z 


n =3 


^ik x Ax+ik y Ay 


y dk x d k y (2.3) 


where Ax - x - x 1 , Ay - y - y', Az = z- z', u(-) denotes the Heaviside step function, and {e p ,„, k Pt „ z , a pn } stand for the 
electric field polarization state vector, longitudinal wave number component, and direct field polarization amplitude 
of the /zth formation bed’s mh plane wave polarization (1 < p < N, 1 < n < 4), respectively. Now introducing 
additional formation beds will induce a modification, via reflection and transmission mechanisms interfering with the 
direct field, to the total observed electric field. We mathematically codify this interference phenomenon by deriving 
the (transmitter layer [ M] and receiver layer [L] -dependent) time-harmonic scattered electric field £ v (r) |fl9l : 


£,( r) = 


(2 n) 2 


+ °° [2 4 

ff (1 - 5l, N ) £ «l„e L ,«e' W + (1 - S ul ) £ «l n e L ,„e^» 

^ Z. n= 1 n =3 


x e 


ik x Ax + ik y Ay dkxdk ( 2 . 4 ) 


where ci s pn is the scattered field polarization amplitude of the nth polarization in layer p, and 5p\_pi denotes the 
Kronecker Delta function. 


2.2. Tilted Layer Modeling 

Admit an /V-1 aver medium where the mth planar interface (m=l,2,...,N - 1) is characterized as follows. First, 
its upward-pointing area normal vector z' m is rotated by polar angle -90° < a' m < 90° relative to the z axisQ and 
azimuth angle 0° < [i' m < 180° relative to the x axis. Second, the mth interface’s “depth” z' m is defined at the Cartesian 
coordinate system’s transverse origin (x,y) = (0,0). See Figure |T] for a schematic illustration of the environment 
geometry’s parametrization. 


2 eo» c, and fJo = l/(eoc 2 ) represent vacuum electric permittivity, vacuum speed of light, and vacuum magnetic permeability, respectively. 
Furthermore, w = 2 nf is the angular temporal radiation frequency, ko = cj/c is the vacuum wave number, and t/o = yfpo/ro is the intrinsic vacuum 
plane wave impedance l‘34jl6il . 

’Albeit as becomes apparent below, this “polar” angle corresponds to rotation in the direction opposite to that ascribed to the spherical coordinate 
system. 
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Figure 1: (Color online) Figure [Ta] shows the original problem with tilted planar interfaces in an TV-layer geological formation 
possessing the EM material tensors {e p ,p. p , & p ). Figure[Tb]shows the transformed, equivalent problem obtained through employing 
special “interface-flattening” media (c.f. Eqn. t |2.7[ l) that coat the underside ({c^ +1 ,ju', +1 ,0’', +1 }) and over-side ({e",/*", <x"j) of the 
mth interface, d represents the thickness of each T.O. slab in meters. For simplicity of illustration, all interfaces here are tilted 
within the xz plane (i.e., interface-tilting azimuth orientation angles {/?'„) = 0°). 


To make the /rath planar interface parallel to the xy plane yet retain its tilted-interface scattering characteristics, as 
the first of two steps we abstractly define, within the two slab regions ([z' m - cl] < z < z' m and z' m < z < [z' m + d]) 
bounding the mth interface, a “coordinate stretching” transformation. Namely this transformation relates Cartesian 
coordinates ( x,y,z ), which parametrize the coordinate mesh of standard “flat space”, to new oblique coordinates 
(x,y,z) that parametrize an imaginary “deformed” space whose coordinate mesh deformation systematically induces 
in turn a well-defined distortion of the EM wave amplitude profile within said slab regions Il54l l55l [62 1 1^] 

x = x, y -y, z = z + a m x + b,„y (2.5) 


where a m = - tan a' m cos [i' m and b m = - tan a' m sin fj' m . Indeed this coordinate transform will cause wave fronts to inter¬ 
act with the mth flattened interface as if said interface were geometrically defined by the equation z, = z' m - a m x - b m y 
rather than z = z' m . As the second step in the interface-flattening procedure, we invoke a “duality” (not to be con¬ 
fused with duality between the Ampere and Faraday Laws) between spatial coordinate transformations and doubly- 
anisotropic EM media which “implement” the effects, of an effectively deformed spatial coordinate mesh (and hence 
effectively deformed spatial metric tensor), on EM waves propagating through flat space (see references deriving this 
“duality” ll471l54H57l ). Following one of two common, equivalent conventions Il47ll54ll leading seamlessly from coor¬ 
dinate transformation to equivalent anisotropic material properties, by defining the Jacobian coordinate transformation 
tensor El: 
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within the region ( z' m - d) < z < z' m one has the interface-flattening material tensors {y' m+l } in place of the original 
formation’s material parameters {f m+ 1 } within layer p = m + 1 (y = e,p,cr). Similarly, within the region z' m < z < 


4 We remark in passing upon a strong similarity between the coordinate transform shown in Eqn. {23} versus the “refractor” and “beam 
shifter” coordinate transforms prescribed elsewhere 16211631 . However, while the beam shifter transform (and the equivalent anisotropic medium 
it describes 1621 ) is perfectly reflectionless due to continuously transitioning the coordinate transformation back to the ambient medium (e.g., free 
space), our coordinate transformation is inherently discontinuous. Indeed, note in Eqn. {23} that the mapping z , which depends on x and y in 
addition to z, can not be made to continuously transition back to the (identity) coordinate transform z(z) = z implicitly present within the ambient 
medium. Alternatively stated, our defined anisotropic coating slabs have the exact same material properties as the beam shifter, but our slabs border 
the ambient medium at planes that, though parallel to each other, are orthogonal relative to the junction planes between the beam shifter and its 
ambient medium. See other references for the importance of the junction surface’s orientation in ensuring a medium perfectly impedance-matched 
to the ambient medium I641I65I . 
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( z' m + d) one has the interface-flattening material tensors { 7 "} in place of the original formation’s material parameters 
{fm} within layer m. How are the interface-flattening material tensors defined though? Quite simply, in fact, and 
this definition holds regardless of the original formation layer’s anisotropy and loss (“T” superscript denotes non- 
Hermitian transpose) (54): 

f m+ 1 = A-m ' fm+ 1 ' A m ; f m = A m ■ f m ■ A m (2.7) 

Note that if the mth interface lacks effective tilt then A m reduces to the identity matrix, which in turn leads to the 
two interface-flattening media reducing to the media of the respective formation layers from which they were derived 
using Eqn. £3= r m+ 1 = 7 „, + i and 7 " = y,„ (as expected!). Now the new material profile, characterized by parallel 
planar interfaces, appears for a simple three-layer, two-interface geometry as: 


71 , (zj + d) < z < °° 

( 2 . 8 ) 

r", z\ < z < (zj + d) 

(2.9) 

72 , (zj - d) < Z < zj 

( 2 . 10 ) 

72, (z ' 2 + d) <Z< (zj - d) 

( 2 . 11 ) 

?2 , z 2 <z<(z 2 + d) 

( 2 . 12 ) 

7 j, ( z 2 -d) <z<z 2 

(2.13) 

73 , -oo < Z <(z 2 - d) 

(2.14) 


with an analogous material profile resultant for N > 3 layers. 

There are two advisories worth mentioning. First, we recommend adaptively (i.e., depending on the transmitter 
and receiver positions) reducing the thickness d of coating layer(s), within which receiver(s) and/or transmitter(s) 
may reside depending on their depth, just enough so that the receivers and transmitters are located once more within 
the formation layers. Why this recommendation? Although pseudo-analytical techniques are available to compute 
fields when the receiver and/or transmitter are located in such layers ( 66 ), the main reason is to eliminate spurious 
discontinuities of the normal ( z in our case) electric and magnetic field components manifest when the source (or, 
as can be anticipated from EM reciprocity, the receiver) traverse a boundary separating a true formation layer and a 
coating slab (47). Second, the thickness d of the coating slabs must also be adjusted to ensure the coating layer just 
beneath the mth interface does not cross over into the coating layer just above the (m + l)st interface. We account for 
these two points within our numerical results below. 

3. Error Analysis 

3.1. Overview 

To briefly recap: The proposed method relies upon insertion of specially-prescribed, doubly-anisotropic material 
slabs above and below each (effectively) tilted original interface to manipulate wave-fronts such that they interact with 
the coated interfaces as if they were tilted. This technique allows one, in principle, to unequivocally and independently 
prescribe the arbitrary, effective polar and azimuthal tilting orientation of each interface. Moreover, in the limit 
of vanishingly small effective polar tilt for some interface (and irrespective of effective azimuth tilt), the material 
properties of the slab just above (below) this interface continuously transitions back to the material properties of the 
ambient medium just above (below) said interface. However, as the material properties of each coating slab are (for 
non-zero effective tilt) not perfectly impedance matched to the respective ambient medium into which it is inserted, 
spurious scattering will result that coherently interferes with the true field scattered from the effectively tilted interface. 

The spurious scattering corrupts the true responses (both co-polarized and cross-polarized ones) arising from 
tilt and hence needs to be quantified if we are to understand the practical limitations of the proposed tilted-layer 
algorithm, which is inherently perturbational in nature with respect to the range of (effective) polar tilt that can 
be reliably modeled. It is also important to understand the error trends not just versus effective polar tilt, but also 
for different geophysical media (anisotropy and loss), transmitter/receiver spacing, sensor positions (relative to the 
interface), transmitter and receiver orientation, and thickness d of the coating slabs. To simplify the geometry and 
admit a closed-form reference solution, we examine only a two-layered medium with an interface effectively tilted 
within the xz plane by a - arj degrees. Both the reference (/7 R / f ) and algorithm/transformed-domain (7/™) field values 


6 


are computed to fourteen Digits of Precision (DoP), with the relative error between their respective field values (either 
real or imaginary part) denoted e (-Log 10 e denotes DoP agreement): Either e = |Re[//™] - Re[ H ^ q ' ]|/1Re[ 1 ] or 
e = |Im[//™] - Im[ffJ^]|/|Im[Hj^]|, where Re[A/ HY/ ] = H' wq and Im [H wq ] = H" q denote the real and imaginary parts 
(resp.) of the ^-oriented magnetic field component observed at the receiver due to the w-directed magnetic dipole 
transmitter H wq ( w, q = x, y, z). Errors below 10 14 were artificially coerced in post-processing to 10 14 . 

In order to compute the reference field solution in this two-layer scenario, we first denote the transmitter-receiving 
spacing L s , transmitter depth z', and receiver depth z - z' + L s in the “transformed” domain with a flat, parallel (to 
the xy plane) interface (residing at z\ = 0) with two coating slabs]-’] In the equivalent domain, involving again a flat 
interface (z\ = 0) but with a rotated sensor, the new transmitter position (x' T ,y' T = 0 ,z' T ) writes as z' T = z! cos a and 
x’ T = z r sin cc, the new receiver position (xr,yr = 0, zt) writes as zt = (z' + L s )cosa and xj = (z! + L J )sina, the 
transmitting dipole orientations are physically rotated by —a degrees in the xz plane, and the receiver antennas are 
also physically rotated by -a degrees in the xz planej^] Additionally, for the latter two material profiles (described, 
and denoted M3 and M4, below) involving anisotropic media, we rotate the anisotropic material tensors by -a de¬ 
grees (isotropic tensors are invariant under rotation, by definition). Throughout this error study, when computing the 
reference and transformed-domain solutions the transmitter radiation frequency is held fixed at f- 100kHz. Further¬ 
more, both the reference and transformed-domain solutions are computed using the same numerical code, albeit with 
effective interface tilt set to zero for the reference solutions. 

We show relative error Log 10 e[//', 9 ] and Lcg 10 e[//" ? ] versus the effective polar tilt angle a < 0 (degrees) and 
azimuth angle (fixed at /?' = 0°). To better illustrate error trends (e.g., linear or quadratic) versus a, we plot this Log- 
scale error on the vertical axis and Logio[a| on the horizontal axis. To show error variations versus material profile 
and transmitter-receiver spacing, in each plot we exhibit relative errors for two different sensor spacings L s (=400mm 
[“SI”] and =lm [“S2”]) and four different conductivity profiles (each having two layers, with e r — /./,• = 1, prior to 
inserting the two coating slabs): 

1. <T\ = lmS/m, <T 2 = 2mS/m ([“Ml”]- Both layers isotropic, 2:1 conductivity contrast) 

2. cri = lmS/m, cr 2 = 20mS/m ([“M2”]- Both layers isotropic, 20:1 conductivity contrast) 

3. cr 1 = lmS/m, cr/a = 5mS/m, cr V 2 = lmS/m, az = 60°, @2 = 0° ([“M3”]- Deviated uniaxial bottom layer) 

4. cr 1 - lmS/m, cr x 2 = 5mS/m, cr y 2 = 2.5mS/m, cr Z 2 = lmS/m, a 2 = Pi = 0° ([“M4”]- Non-deviated biaxial 
bottom layer) 

where the anisotropic tensor components (<t /,2 , cr v 2 , cr x 2 , cr v 2 , cr z 2 1 and tensor spherical (polar and azimuth) deviation 
angles {cx 2 ,Pi) define the principal conductivity values and orientation (resp.) of the conductivity tensor relative 
to the cartesian axes (discussed in Section [4] and |j2]). The curve labeled “SXMY” refers to the sensor spacing 
denoted above in square brackets by [“.S’ A - ”] (X=l,2) and the material scenario denoted above in square brackets by 
[“MY”] (7=1,2,3,4). Fixing the flattened interface’s depth at z\ = 0, for two coating slab thicknesses (cl =2mm 
[“dl”] and d = 200mm [“d2”]) we examine three mid-point sensor depths in the transformed domain: D = 2m 
(transmitter and receiver in top layer [“Ol”]), D = 0m (transmitter in bottom layer, receiver in top layer [“02”]), and 
D = -2m (transmitter and receiver in bottom layer [“03”]). For the real and imaginary part of each field component 
we examine six sensor-location/t/ permutations, which are denoted on each page with plot labels (“Scenario:01dl”, 
etc.). The transmitter depth (z r ) and receiver depth (zj are computed using the sensor spacing L s and mid-point depth 
D: z = D + L s /2 and z' = D - LJ 2. 

A final remark before proceeding with error analysis, concerning our choices of examined d: One could in princi¬ 
ple make the slabs arbitrarily thin or thick, which we have not tried (we only examined, in this error study, d- 0.2m and 
of=2mm). The algorithm’s present design does not dictate a specific “optimal” value of d, unfortunately. What we can 
say however is that making d comparable to the local wavelength, on either side of the interface, is highly undesirable 
for at least two reasons. First, since the coating slabs are in fact reflective they can confine spurious guided-wave 
modes that may significantly corrupt computed sensor responses. Second, as indicated in Section 2.2, the thicker d 
is one must adaptively reduce the thicknesses of slab(s) intersecting each other, as well as slab(s) containing trans¬ 
mitters or receivers, to eliminate artificial discontinuities in the normal field components. This adaptive thickness 


5 Note: For all results in this paper, we assume a vertically-oriented sensor in the transformed domain. 

6 Equivalently, the observed fields at the receiver location are now observed relative to a cartesian coordinate system rotated by —a degrees. 
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reduction, which becomes increasingly frequent for finer spatial sampling (i.e., versus sensor depth) of the sensor re¬ 
sponse profiles, introduces yet another level of arbitrariness into the algorithm. Namely, the desirable minimum space 
maintained between the receiver and coating layer/ambient medium interfaces (and likewise for the transmitter), as 
well as between any two ambient medium/coating slab junctions. Although both examined values of d in this error 
study would not necessitate their adaptive thinning, due to the three specifically examined sensor locations, the thicker 
d =0.2m slabs would require adaptive thinning in the Section [4] results due to the sensor’s depth being varied (0.1m 
sampling period) throughout the studied three-layer formation profiles. 

For the cross-polarized field plots, we only show errors for H xz and H zx since the other cross-polarized field com¬ 
ponents (H xy , H vx , H yz , H zy ) had zero magnitude to within numerical noise. Elaborating: Our chosen threshold for 
numerical noise is that either 1(T 12 and/or | H™\ < 10 12 (-12 on Logio scale), which is based on the adaptive 

integration tolerance (1.2 x 10 I4 ) and an educated guess (10 2 ) of the maximum magnitude of field components that 
would experience cancellation during evaluation of the oscillatory Fourier double-integral (Eqns. In 

passing, we mention that to suppress numerical cancellation-based noise due to using finite-precision arithmetic, we 
“symmetrically” integrate. That is to say, we integrate along the integration contour sub-sections (a, b ) and (- b , -a), 
add these results, then integrate along contour sub-sections (b, c) and (-c, add these results and update the accu¬ 
mulated contour integral result, and so on. 


3.2. Results and Discussion 

Figures[2pT|display errors (respectively) for the following field components: Re[77 vv ], Im ['H xx ], Re[77 vv ], Im[77 vv ], 
Re [77-,], Im[77,,], Re[77/.-], Im[77 A ,], Re[77, A ], and Im[77-J. Let A f denote the rate of change of error (Log 10 e) versus 
Log HI |«|: A slope of two (one) indicates quadratic (linear) error variation versus effective tilt. We remark that in some 
plots, one or more material/tool-spacing scenario curves show strange “dips” in the error levels (Figs. |2a||2b| [5a] |6d] 
m and[8a]). Given how small the error dips are (typically < 1 DoP), we ascribe the dips to a combination of machine- 
dependent computation and problem geometry-dependent numerical cancellation (effective tilt, material profile, d, 
and sensor characteristics). There are also some “kinks” in the error behavior, at very small tilt, that can be observed 
in one or more material/tool-spacing scenario curves within each of Figures 2cpf 3c| 3e |4c]|4fj 5aJ 5c 5e 6cfl6f 
[7c] and[7e] Despite these two sporadically-occurring characteristics in the error curves (typically only 1-2 curves in 
each stated figure has these properties), the overall error trends are well preserved, and it is this we summarize in the 
following observations: 


1. A two order of magnitude increase in d (from 2mm [“dl” plots] to 200mm [“d2” plots]) produces low error 
variation (typically 0-1 DoP error increase). 

2. Error is typically 1-3 DoP higher when the transmitter and receiver are in different layers (“02” plots), versus 
when they are in the same layer (“Ol” and “03” plots). By contrast, the error is approximately equal if the 
transmitter and receiver are both either above or below the interface. 

3. Transmitter-receiver spacing (“SI” vs. “S2” curves) has little effect on error levels (typically 0-1 DoP differ¬ 
ence). 

4. The M3 scenario curves typically show greatest error (versus Ml, M2, and M4 curves) in the co-polarized 
results, but the lowest error in the cross-polarized plots. There does not appear to be any obvious, systematic 
trend in error variation between the Ml, M2, and M4 cases across the studied sensor/environment parameter 
permutations. 

5. The Ml, M2, and M4 curves show quadratic error variation in the co-polarized results (their cross-polarized 
errors are intolerably high), while the M3 curves predominantly show instead linear error variation for both co- 
polarized and cross-polarized results. For Figures |3c||3d| [5c|5d| and [7c||7d1 the M3 curves, interestingly, show 
quadratic error too however. 

6. The cross-pol response errors not only are much higher than their co-pol response error counterparts, but the 
errors vary more versus d, L s . and sensor position D too. 


7 Although rigorous justification for the threshold 10 17 is lacking, the conclusion of negligible {!l xy . //,,, H zy ) is physically reasonable. 

Indeed one can verify (using Eqns. (2.6) - ]2i7] ) that T.O. media, used to effect modeling of xz plane tilting, only perturbs the ambient medium’s 
response in the xz plane (but not v direction). Hence one can reason the T.O. media should not induce spurious scattering leading to non-zero 
j H^, H yx , //,-, H z y) responses (i.e., the cross-pol responses involving ay-directed transmitter or receiver) if they were absent without the T.O. slabs. 
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Figure 2: Relative error in computing H' xx =Rs\H xx ], 
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Figure 3: Relative error in computing H" x =\m\H xx \. 
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Figure 4: Relative error in computing H' yy =Re[H yy ]. 
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Figure 5: Relative error in computing H yy =\m\H yy \. 
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Figure 6: Relative error in computing //l_=Re[//--]. 
13 

































H ]3 ut 6on L H]3 Ut 6on L, H]3 6o-| 


Scenario:01d1 


Scenario:01d2 



Logjal 

(a) 


I 

UJ 


O) 

o 



Scenario:02d1 


Scenario:02d2 



(c) 


(d) 


Scenario:03d1 


Scenario:03d2 




Log 10 |a| 


-1 -0.5 0 

Log 10 M 


(e) 


(f) 


Figure 7: Relative error in computing H”=\m\H zz \. 
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Figure 8: Relative error in computing H' xz =Rs\H xz \. 
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Figure 9: Relative error in computing 
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Figure 10: Relative error in computing H' zx =Re[H zx ], 
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Figure 11: Relative error in computing H" x =\m\H zx \. 
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4. Application to Triaxial Induction Sensor Responses 


Now we perform case studies, to illustrate the algorithm’s flexibility in modeling media of diverse anisotropy and 
loss, involving twenty eight variations of a three-layered medium (seven tilt orientations, four central bed conductivity 
tensors), where the two interfaces exhibit (effective) relative tilting and are each flattened by two coating slabs d- 2mm 
thick (i.e., one coating slab immediately above, and one coating slab immediately below, each flattened interface). 
Figure [TJJdepicts the well-logging scenario simulated: A vertically-oriented triaxial induction tool {7][8l, that operates 
at f- 100kHz, is transverse-centered at (x, y) = (0,0), possesses three (co-located) electrically small loop antenna 
transmitters (modeled as unit-amplitude Hertzian [equivalent] magnetic current dipoles) directed along x ( M T X ), y 
( M R ), and z (Ml), and possesses three (co-located) loop antenna receivers \M R , M K , M R ) situated L s = 1.016m (forty 
inches) above the transmitters]^] 

All four material scenarios share the following common geological formation features prior to inserting the 
interface-flattening slabs: Top interface depth z\ = 2m and bottom interface depth z ' 2 = -2mj^all three layers possess 
isotropic relative dielectric constant (i.e., excluding conductivity) e r \ = e r2 = e ,-3 = 1 and isotropic relative magnetic 
permeability /j r \ = /i r2 = fJ,z = 1, layer one (top layer) possesses isotropic electric conductivity <X] = 50mS/m, layer 
three (bottom layer) possesses isotropic electric conductivity cr 3 = 20mS/m, the polar tilt angles of the two interfaces 
are equal to a' 2 = —a[ — a' > 0, and the azimuth tilting orientation of the two interfaces are equal to /?' - — /f. 

In other words, the relative tilting between the two interfaces is 2a' degrees while we choose to tilt both interfaces 
with identical azimuth orientation. This choice is made to isolate and better understand effects of relative polar tilting, 
without the added, confounding effect of relative azimuth deviation between the interfaces. 

In the Figures below, the labels Tl, T2, etc. in the legend represent induction log signature curves corresponding 
to different interface tilting scenarios, namely: {a' = 0°,/T = 0°) (Tl, solid black curve), {a' = 1°,/?' = 0°) (T2, solid 
blue), {a' = \°,p = 45°) (T3, solid green), {a' = 1 °,/3' = 90°) (T4, solid red), {o' = 3°,/T = 0°) (T5, dash-dot blue), 
(a' = 3°,/i' = 45°} (T6, dash-dot green), (a' = 3°,/?' = 90°) (T7, dash-dot red). We chose the azimuth tilt orientations 
/?' = 0° and (j = 90° specifically due to our having examined already, in the previous section, the error of co-polarized 
fields components when they are oriented either within or orthogonal to the plane of interface tilting. Moreover for a 
given a', the /?' = 45° curves exhibit results intermediate to the // = 0° and jj = 90° results, suggesting confidence in 
the /?' = 45° results too. 

(Material) Scenario 1: The highly resistive central layer (mimicking a hydrocarbon-bearing reservoir) possesses 
isotropic electric conductivity <T 2 = 5mS/m. Scenario 2: The central layer’s conductivity is characterized by a non- 
deviated uniaxial conductivity tensor d- 2 = o 7 !2 (xx + yy) + cr v2 zz with <t/ i2 =5mS/m and cr v2 =lmS/m. Scenario 3: 
Same as Scenario 2, except the symmetric, non-diagonal conductivity tensor writes as [2} 


tr 2 


( T xx2 

&xy2 

°~xz2 

0~xy2 

0-yy 2 

0~yz2 

_Vx Z 2 

Cyz2 

0~zz2 


o- X x2 = tr h 2 + (cr v 2 - cr ,, 2 )(sin a 2 cos p 2 f 
crxyl = ( cr v2 - cr h 2)(sma 2 ) 2 sin/? 2 cos p 2 
cTxzi = (o \2 - CT; a) sin a 2 cos a 2 cos [h 
CTyy2 = cr h 2 + (cr v 2 - (r h 2 )(sina 2 sin / 3 2 ) 2 
o- yz 2 = ( cr V 2 - cr h2 ) sin a 2 cos a 2 sinj 6 2 
o~zz 2 = cr v2 - (cr v2 - cr A2 )(sin a 2 ) 2 


(4.1) 

(4.2) 

(4.3) 

(4.4) 

(4.5) 

(4.6) 

(4.7) 


with tensor dip (a 2 ) and strike (P 2 ) angles equal to a-i = 30° and p 2 = 0° (compared to Scenario 2, where 01,2 = Pi = 
0°)JTo] Scenario 4: The central layer has full (albeit non-deviated) biaxial anisotropy characterized by the conductivity 
tensor &2 - cr^xx + cr v2 yy + cr z2 zz, where {cr v2 = 5, cr y2 = 20, cr z2 = l}[mS/m]. 


8 Note 1: “T” superscript here denotes “transmitter”, not non-Hermitian transpose. 

9 For the equivalent problem with flattened interfaces, this choice of z[ and z' 2 results in the coating layer just above the bottom interface, and 
coating layer just below the top interface, having their respective coating layer/central formation layer junctions spaced 3.996m apart. 

10 Note: These conductivity tensor dip and strike angles, completely unrelated to the interface tilting angles describe the material 

tensor’s polar and azimuthal tilting (resp.) but follow a different convention (3). 
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Now we discuss the co-polarized results; note that as the cross-polarized results cannot be reliably modeled nearly 
as well as the co-pol fields (see previous section), we omit them. Observing Figure [13] 

1. The real part of all three co-polarized measurements has no visibly noticeable sensitivity to interface tilting, 
even when the sensor is near bedding junctions. The imaginary part of these measurements does, by contrast, 
exhibit tilting sensitivity. 

2. The imaginary parts’ sensitivity to tilting depends on the sensor position, with sensitivity being largest when the 
sensor is near the interfaces (+2m). This is expected, since the conductive formation exponentially attenuates 
fields scattered and propagating away from the interfaces. 

3. The 4m bed thickness and sensor frequency (100kHz) preclude observation of inter-junction coupling effects 
due to wave “multi-bounce” within the slab layers ll34l rCh. 2]. 

4. At a fixed polar interface tilting (T2, T3, and T4 [±1°] versus T5, T6, and T7 [+3°]), Im[//-,<,/] shows no visible 
sensitivity to the interface’s azimuth orientation. This lack of azimuthal sensitivity is quite sensible, given the 
azimuthal symmetry of this z-transmit/z-receive measurement a. 

5. In contrast to Im[fly], Im[fly] and Im[//,,y ] do show azimuthal sensitivity. For fixed polar tilt, both measure¬ 
ments show greatest excursion (i.e., away from the black, zero-tilt curve) when the interfaces’ common azimuth 
tilt orientation is aligned with the transmitter and receiver orientation (0° [blue curves] and 90° [red curves] 
for Im[fly] and Im[fl v y], respectively). On the other hand, minimal excursion occurs when the azimuth tilt 
orientation is orthogonal to the transmitter and receiver orientation (90° and 0° for Im [H X ' X '] and Im[fl v y], 
respectively). 

6. Tilting, irrespective of the polar orientation’s sign (e.g., 1° versus -1° polar tilt), results in the responses 
Im[fly] and Im[fl v y ] having downward excursions. Indeed for a fixed azimuth tilt (curve color), observe 
the sensor response near the two (effectively) oppositely-tilted interfaces for both the solid and dash-dot curves. 

7. In contrast to Im[//,v] and Im[fl v y ] exhibiting downward excursions, Im[fly ] always shows an upward 
excursion irrespective of polar tilt sign. 

The above observations also visibly manifest for the three remaining cases involving anisotropic media. We find, for 
the small tilting range explored here at least (zero to six degrees of relative tilt), that the anisotropy primarily serves 
to alter “baseline” (zero-tilt) sensor responses which are then perturbed by the effect of tilt. 
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Figure 12: Original geometry (Fig. 12a i and transformed, approximately equivalent geometry (Fig. 12b i employed below. For 
clarity of illustration, the layers are shown tilted within the xz plane (j3' =0°). 
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Figure 13: Co-polarized, complex-valued received magnetic fields: Material Scenario 1. 
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Figure 14: Co-polarized, complex-valued received magnetic fields: Material Scenario 2. 
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Figure 15: Co-polarized, complex-valued received magnetic fields: Material Scenario 3. 
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Figure 16: Co-polarized, complex-valued received magnetic fields: Material Scenario 4. 
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5. Conclusion 


We have introduced and profiled (both quantitatively and qualitatively) a methodology, which augments robust 
full-wave numerical pseudo-analytical algorithms admitting formation layers of general anisotropy and loss, to incor¬ 
porate the effects of planar interface tilting. Our proposed methodology, directed at extending the range of applica¬ 
bility of such eigenfunction expansion methods by relaxing the traditional constraint of parallelism between layers, 
consists of first defining a spatial coordinate transformation within a thin planar region surrounding each interface 
to be “flattened”. Such coordinate transformation locally distorts the EM field, effectively altering the local angle 
of incidence between EM waves and flat interfaces so as to mimic the presence of interfaces possessing effective, 
independently-defined tilt orientations. During the second stage of flattening the interfaces, the coordinate transfor¬ 
mation is incorporated into the EM material properties of said planar regions via application of T.O. principles, which 
exploits the well-known “duality” between spatial coordinate transformations and equivalent material properties “im¬ 
plementing” these coordinate transformations in flat space. 

As the proposed methodology is not limited by the loss and anisotropy properties of any layer, some combina¬ 
tion of beds such as non-deviated sand-shale laminates, cross-bedded clean-sand depositions, formations fractured 
by invasive drilling processes, and simpler isotropic conductive beds can all be included along with flexibly-defined 
(effective) bed junction tilting with respect to polar deviation magnitude as well as azimuth orientation. Exhibiting 
representative examples of the new methodology’s said flexibilities, we applied it to demonstrating the expected qual¬ 
itative properties of multi-component induction tool responses when planar bedding deviation is present. This being 
said, one should not confuse modeling flexibility and manifestation of expected qualitative trends with quantitative 
accuracy. Indeed, although T.O. theory informs us that said interface-flattening media facilitate numerical predic¬ 
tion of tilted-interface effects via employing specially-designed coating slabs, the necessary orientation of the slabs’ 
truncation surfaces (i.e., planes parallel to the xy plane) predicts that spurious scattering will arise. Particularly, our 
numerical error analysis demonstrated that artificial scattering typically scales quadratically with effective interface 
tilt. We found that one could nonetheless safely model tilting effects so long as the magnitude of each interface’s polar 
deviation is kept small, which guided the choice of explored interface tilt ranges when examining qualitative trends of 
sensor responses produced by induction instruments within tilted geophysical layers. 
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